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' Summary. The metabolic networks are very well characterized for a large set of 

organisms, a unique case in within thelarge-scale biological networks. For this reason 
they provide a a very interesting framework for the construction of analytically 
tractable statistical mechanics models. In this paper we introduce a solvable model 
for the distribution of fluxes in the metabolic network. We show that the effect of 
the topology on the distribution of fluxes is to allow for large fluctuations of their 
^ I values, a fact that should have implications on the robustness of the system. 

, 1 Introduction 

Dynamical models on networks have attracted a large interest because of the 
^ \ non-trivial effects of network structure [1, 2, 3, 4] on the dynamics defined 

on them [5] . Important examples of the dynamics on networks with relevant 
QQ \ applications are the Ising model [6, 7, 8], the spreading of a disease [9] and the 

' synchronization models [10, 11]. In this paper introduce a solvable model for 

jy-^ , the distribution of fluxes in the metabolic network. While motivations come 

' from the study of the metabolic network, the problem is quite general and 

, can be applied to supply networks and to many other linear problems [12] of 

' constraint satisfaction on continuous variables on a network. 

Metabolic networks describe the stoichiometric relations between sub- 
. J ■ strates in biochemical reactions inside the cell. They have been mapped [13] 

r> \ for a large number of organisms in the three different domains of life (archea, 

' bacteria and eukaryotes). They provide the biomass needed for cell duplica- 

tion, and the rate of biomass production (growth rate) can be identified with 
a fitness of the cell. The structure of the metabolic network can be repre- 
sented as a factor graph with nodes that are chemical reactions and function 
nodes that are chemical metabolites. The projection of the network on the 
metabolites has a power-law degree distribution and a hierarchical structure 
[14, 15, 16]. To each factor node, which describes a chemical reactions, it is 
associated an enzyme which itself is produced by a regulated gene network. 
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Important aspect of the functioning of these very complex systems include 
dynamical considerations. Flux-balance-analysis [17, 18, 19] make a major 
semplification in the problem. In fact it considers only the steady state of the 
dynamics and includes all the dynamical terms inside the definition of the flux 
of a reaction. For this reason it was able to predict with sufficient accuracy 
the fluxes of the reactions in the graph for a given environmentand it con- 
situte a real break-through in the field. Special interest has been addressed 
to the perturbation of the distribution of the fluxes after knockout of a gene 
or in different environments [20, 21]. The problem of identifying the flux dis- 
tribution in Escherichia coli was studied experimentally [22] and by means 
of Flux-Balance- Analysis [23] . A fat tail in their distribution with different 
power-law exponents a <2 was found. 

Metabolic networks provide a very interesting framework for the construc- 
tion of analytically tractable models using tools of statistical mechanics of 
disordered systems. In this paper we will discuss the impact of the network 
structure (degree distributions) on the steady state distribution distribution 
of the fluxes. Wc shall consider random networks with the same degree dis- 
tribution as the real ones i.e. networks in the the hidden-variable ensemble 
[24, 25, 26] with same expected degree distribution as the metabolic factor 
graphs. Formally the problem is resolved with replica calculations on diluted 
networks [7] extended to the case of continuous variables. Due the simplicity 
of the Hamiltonian the problem is solved with an expansion of the order pa- 
rameter in terms of Gaussians. The problem shares some similarity with other 
problems in statistical mechanics of disordered systems [28, 29]. In a recent 
paper [30] a similar problem was considered in the framework of a different 
model where the steady state of the fluxes is not a priori considered and the 
positive fluxes don't have any upper limit. 



The metabolic; network has a bow tie structure [16] : the metabolites are 
divided into: (i) input metabolites which are provided by the environment, (ii) 
the output metabolites which provide the biomass and (Hi) the intermediate 
metabolites. The stochiometric matrix is given by where ii — 1, . . . , M 

indicates the metabolite and i = I, . . . , N the reaction and the sign of ^^^^ 
indicates if the metabolite /i is an input or output metabolite of the reaction 
i. As in the Flux-Balance- Analysis method we assume that each intermediate 
metabolite has a concentration which is consumed/produced by a reaction 
i at a rate /j. At steady state, we have 



2 The model 




(1) 



where fi is the flux of the metabolic reaction i. For the metabolites present 
in the environment and the metabolites giving rise to the biomass production 



Viable flux distribution in metabolic networks 



3 



we can fix the incoming flux given by 

i 

Th.e fluxes in Eqs.(l-2) can vary inside a fixed volume J?. We assume for sim- 
plicity that this volume is an hypercube fl = [0, 2L]^. Changing the variables 
/j in the variables Si = fi — L and the equations that the fluxes Sj must satisfy 
are given by 

N 

XI ^i^'i^i = Quioi: iJ,= l,...M. (3) 

i=l 

where = a'^ — L Cm,*- "^he volume of solutions V, given the constraints 
(3), is proportional to the quantity 

''0 ■' i=l i j 

where we have used the heterogeneous spherical constraints 

and integrated over L' in the interval [0, L] in order to allow analytical treat- 
ment of the problem. 



3 Replica method 

We assume that the support of our stochiometric matrix is a random network 
with given degree distribution, i.e. a realization of the random hidden-variable 
model [24, 25, 26]. In particular we fix the expected degree distribution of the 
nodes of the factor graphs to be qi for the reaction node i = 1,. . .N and for 
the metabolite nodes ^ = 1, . . . , M and we assume that the matrix elements 
^n^i are distributed following 

= [^(^-^ - + ^^^-^ + ^{'-^) 

where i5() indicates the Kronecker delta. Note that in (6) we have assumed 
that the elements of the stochiometric matrix have values 0, ±1 with a random 
sign and that the variables qi qf^ are nothing else than the hidden-variables 
associated with metabolite n of the reaction i of the hidden- variable network 
ensemble [24, 25, 26]. 

In order to evaluate the steady state distribution of the fluxes in a typical 
network realization we replicate the realizations of the and we compute 
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(log(V)) over the different network realizations. To calculate this average we 
use the replica trick S = (log(Z)) = lini„_>o J""^ . The averaged unormal- 
ized volume of solutions < > can be expressed as 

f-L 



<V">= dL' J II duj"' J II ds,,a J n exp 



exp 



•(7) 



Using the techniques coming from the field of diluted systems, we introduce 
the order parameters [31, 7] 



(8) 



getting for the volume 



V'j^ J Vc{X) J Vc{X) J Vc{s) J Vc{s)exp[nNE{c{X),c{X),c{s),c{s))] 
with 

nr = y" dXic{X)c{X)+ j dsic{s)c{s) - j dX j dsc(A)c(s)(l - cos(A- s)) + 
+ E j dX exp[-i5^ E ~ i<lA>^)\ - « E ^"^'"^ 



+ 



(5^ 



;^E^°S j dseyiV'[-iqic{s)+i^qiU)asl]. 



The saddle point equations for evaluating S are given by 

ic{X) = I dsc{s) {1 — cos{X- s)) 



ic{s) = j dXc{X){l — cos(A- s)) 



c(A) 
c{s) 



{qi)N 



1 



{qi)N 
1 



E' 
E' 

i 

E^> 



exp [-ig^' J^g - iq^iC(X)] 

lUadK^^pHa^EaK-^Q'-ciX') 

exp [-iqic{s) + i J2a ^a^a] 



/rfs'cxp -iqic{s') +iJ2a'^ais')l 

/ dssl exp[-iqic{s) + i J2a Qi^^sl] 

{qi)N ^ / ds' exp[-iqic{s) + i J2a H^asl] 



(9) 
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We assume that the solution of the saddle point equation is replica symmetric, 
i.e. the distribution of the variables z"" = X°',s'^ conditioned to a vector field 
X are identically equal distributed, 



/n 
dxP{x)'[[4>{z''\x) 
a=l 



(10) 



where (f>(z\x) arc distribution functions of z and P{x) is a probability dis- 
tribution of the vector field x. For the function (/'(z|a;) the exponential form 
is usually assumed in Ising models. In our continuous variable case for our 
quadratic problem, we assume instead that (j){z\x) has a Gaussian form. This 
assumption could be in general considered as an approximate solution of the 
equations (9). Explicitly we assume that the functions c(A) and c(s) can be 
expressed as the following, 



:hsS„ 



c(A) = / dmxdhxP{hx,mx)Y[exp 

^ a 

c{s) = / dmgdhsP{hs,ms)Y[^^P 

(jJa = 

from which we derive for c(s) and c(A) 

c(s) = j dmxdhxP{hx,mx)^exp 

c(A) = —i il — J dmsdhsP{hs, rUg) exp 



1 , ,2 1 m? "I 



C0s[TOAAa] 



1 nig 



cosh[mss"] 



27r 



(11) 



2hx'\ 
2hs 



COSh[TOASa]/ft.Aj 

coslm.X" /hs\ I .(12) 



The saddle point equations (9), taking into account the expression for the 
order parameters (11) (12) closes and can be written as recursive equation for 
P{hx, nix) and P{hs,ms), i.e. 



P{hx,mx) 



At fe 
k 



WdKdm^g^PiK.ml) 



i=l "-sj ^ 



1 

2^ 



2a;g, I^^^J m,-^(-l)-"'^ 



. ni\ 



(13) 
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Finally S can be calculated at the saddle point as 



S = 



dhsdmsdh\dm\P{hs, m,s)P{hx,mx) 
rrismx 



2{hx + l/hs) hs + 1/hx 



+ logcosli 

Xhshx + 1 



(14) 



^1 k 



2k+2 



+ 



4 Population Dynamics 

We solved the equations (13) by a population-dynamical algorithm. We rep- 
resent the effective field distributions {hs, rus) {hx, mx) by a large population 
of P ^ 1 fields. Running the algorithm the population is first initialized ran- 
domly and then equations (13) are used to iteratively replace the fields inside 
the population until convergence is reached. Instead of fixing u) we introduce 
a further variable A fixing the value of the average flux in the network. The 
action of the algorithm is summarized in the following pseudo code 

algorithm PopDyn{{hl,ml, h^, ml,..., hf, mf }; 

{h\, m\, hi, ml,..., h^, }, u) 

begin 

do 



• choose a metabolite io with probability qiP{qi); 

• draw d from a Poisson distribution (e~^*gf /fc!) 

• select d indexes ii, . . .id € {1, . . . M} 

• draw a d-dimensional vector n = {rij} of random numbers Uj = 0, 1 
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1=1 "A 



,*0 

's 




L2 




L2 + 



1 fei° + (mi°)^ 



(15) 



• choose a random reaction /xq with probability QfiPiqn) 

• draw (Z form a Poisson distribution {e~'^i^q^/k\) 

• select d indexes ii, . . .14 S {1, • • • M} 

• draw a d-dimensional vector n = {rii} of random numbers Ui = 0, 1 



while (not converged) 

return 

end 



We run the population dynamics algorithm and we measure the distribu- 
tion of the average fluxes nis/hs, the distribution of the fields hg for different 
values of A. Wo consider as the iinderline network a network with the real 
degree distribution of the metabolic factor graph of Saccharomyces cerevisiae 
and on a network with the same number of metabolites and reactions as the 
real Saccahromyces cerevisiae network but with a fixed connectivity for each 
metabolite and reaction node. We consider a population of P = 3000 pair 
of fields (/isjTOs). A random fraction of 5% of the nodes is chosen as an in- 
put/output metabolite. The values of g'^ are chosen randomly depending if the 
metabolite /x is an intermediate metabolite or an input/output metabolite. 

For the input/output metabolites we assume that g^^ is a random number 
uniformly distributed in the interval [—100/1, lOOyl] mimicking high rate of 
production/consumption. For intermediate metabolites we choose g^ with a 
uniform distribution in the interval [— yl,yl]. 





d 



(16) 
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The distribution of nig/hg as a function of A are plotted in figure 1(a) for 
the real metabolic network of Saccharomyces cerevisiae and in figure 1 (b) for 
the random network with two delta function degree distribution {P{qi) = (qi), 
P((7^) = {qi)N/M). We observe that the average fluxes distribution nis/hs 
in Saccharomyces cerevisiae for low A has a fatter tail for the real degree 
distribution than for the two delta peak degree distribution. 

On the other side the distribution of hs is very different in the real and 
in the random case (see figure l(c)-(d)). In particular for the real metabolic 
network degree distribution is broader allowing with higher probability for 
smaller value of hs than in the case of a two delta peak degree distribution. 
Therefore we have shown that the real topology of the networks has as a major 
conseauence in allowine lareer fluctuations of the fluxes in the network. 




m/h h 



Fig. 1. Distribution of the fields m,s/hs and of the fields hs for a the real degree 
distribution of the metabolic network of Saccharomyces cerevisiae (graphs (a) and 
(c))and for a graph with the same number of metabolites and reactions and the 
same number of nodes that the real metabolic network of Saccharomyces cerevisiae 
but with two delta peaks for the degree distribution (graphs (b) and (d)). 



5 Conclusions 



In this paper we have proposed a statistical mechanics approach for the study 
of flux-balance- analysis in a particular ensemble of metabolic networks. We 
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have studied the impact of the topology of the networks on the distribution 
of the fluxes. Wc observe that the role of the real topological structure is to 
allow for larger variation of the fluxes, a fact that should have implications 
for the robustness of the system. In particular we found that the topology of 
real metabolites has an impact on the fat tail of the nis/hs distribution and 
on the small h field of the network, Further work is under consideration for 
the implementation of a message-passing algorithm able to predict the fluxes 
taking into account the full complexity of the real metabolic network. 
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